function [MaxSR, weights, iSR, RMS, sigmae, maxf]=lettau_rp(K, X, F, eigval)
%This function implements the RP-PCA methodology in Lettau and Pelger (2020a,b)
% It also performs the Onatski test for the number of latent factors
% Input:
% K= number of factors taken into consideration
% X = T x N matrix of portfolio excess returns
% F = T x K matrix of latent factors
% rmax = max number of factors to be tested in the Onatski algorithm for
% the number of factors
% eigval = vector of ordered eigenvalues to perform Onatski algoritmh
% output:
% MaxSR = Max Sharpe ratio obtained by combining the latent factors
% weights = weights attacched to the factors to get MaxSR (This differs
% form iSR as it is normalized by variance and take into account
% correlation between factors (if factors are correlated)
% iSR =individual Sharpe ratio to each factor 
% RMS = square root of the average pricing errors of the TS regressions
% sigmae = unexplained variance implied by the model
% maxf  = optimal number of factors suggested by the Onatski criterion
%% Compute statistics for Lettau and Pelger (2020ab)
T=size(X,1); %time obs
N=size(X,2); %cs obs

for k=1:K    %normalize  factors to have positive average
    if (mean(F(:,k))<0)
        F(:,k)=(-1)*F(:,k);
    end
end 
Bmat=zeros(N,size(F,2));
emat=zeros(T,N);
alpha=zeros(N,1);
xmat=[ones(T,1) F];
for n=1:N
    tempout=inv(xmat'*xmat)*xmat'*X(:,n);
    atemp=tempout(1);
    btemp=tempout(2:end);
    etemp=X(:,n)-(atemp*ones(T,1)+F*btemp);
    alpha(n)=atemp;
    Bmat(n,:)=btemp';
    emat(:,n)=etemp;
end
%% Evaluate performance
RMS=sqrt((alpha'*alpha)/N);
sigmae=sum(diag(cov(emat)))/sum(diag(cov(X)));
sigmaF=cov(F);
muF=mean(F)';
weights=inv(sigmaF)*muF;
Fcomb=F*weights;
MaxSR=mean(Fcomb)/std(Fcomb);
midiovar=sum(diag(cov(emat)))/N;
iSR = muF./diag(sqrt(sigmaF));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Onatski (2010) test
rmax=K;
r=0;
rmaxtemp=rmax;
iter=0;
j=rmaxtemp+1;
while (rmax-r)>0
yo=eigval(j:j+4);
xo=[(j-1)^(2/3) (j)^(2/3) (j+1)^(2/3) (j+2)^(2/3) (j+3)^(2/3)];
xo=[ones(5,1) xo'];
beta=inv(xo'*xo)*xo'*yo;
th=2*abs(beta(2));  
i=j-1;
diff = eigval(1:i,:)-eigval(2:i+1,:);
temp = find((diff-th)>=0);
r=max(temp);
 if isempty(r) 
     r=0;
 end 
if size(diff,1)==r
rmax=rmaxtemp;
end
iter=iter+1;
j=r+1; 
rmaxtemp=rmaxtemp-1;
end
maxf=rmax;
   




